Numerical modeling of a Global Navigation Satellite System in a 
general relativistic framework 

Pacome Delva^ *, Uros Kostic'', Andrej Cadez*' 

"ESA Advanced Concepts Team, ESTEC, DG-PI, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands 
^Department of Physics, University of Ljubljana, Jadranska 19, 1000 Ljubljana, Slovenia 



Abstract 

In this article we model a Global Navigation Satellite System (GNSS) in a Schwarzschild space- 
time, as a first approximation of the relativistic geometry around the Earth. The closed time-like 
and scattering Ught-Uke geodesies are obtained analytically, describing respectively trajectories 
of satellites and electromagnetic signals. We implement an algorithm to calculate Schwarzschild 
coordinates of a GNSS user who receives proper times sent by four satellites, knowing their 
orbital parameters; the inverse procedure is implemented to check for consistency. The con- 
stellation of satellites therefore reaUzes a geocentric inertial reference system with no a priori 
realization of a terrestrial reference frame. We show that the calculation is very fast and could be 
implemented in a real GNSS, as an alternative to usual post-Newtonian corrections. Effects of 
non-gravitational perturbations on positioning errors are assessed, and methods to reduce them 
are sketched. In particular, inter-links between satellites could greatly enhance stability and ac- 
curacy of the positioning system. 

Keywords: General Relativity, Schwarzschild Space-time, Relativistic Positioning System, 
GNSS 



1. Introduction 

The classical concept of positioning system for a Global Navigation Satellite System (GNSS) 
would work ideally if aU satellites and the receiver were at rest in an inertial reference frame. But 
at the level of precision needed by a GNSS, one has to consider curvature and relativistic inertial 
effects of space-time, which are far from being negligible. There are two very different ways 
of including relativity in a positioning system: one way is to keep the newtonian conception of 
absolute time and space, and add a number of post-newtonian corrections depending on the de- 
sired accuracy; another way is to use a relativistic positioning system. This is a complete change 
of paradigm, as the constellation of satelUtes is described in a general relativistic framework. 
This new scheme for positioning potentially allows the definition of a very stable and accurate 
primary reference system. 
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For a detailed review of the post-newtonian scheme of a GNSS one can read Ashby ( 2003| l 
and ,Pascual-Sanchez (2007 1. In this scheme, the two main corrections come from gravitational 
frequency shift between the clocks - due to the local position invariance principle - and from 
the Doppler shift of the second order - due to relative motion of satellites and users. A "GPS 
coordinate time" is defined as the time of a clock at rest on the geoid. To add the post-newtonian 
corrections, it has to be related to the time measured by a clock in an inertial frame at spatial in- 
finity. Then one has to perform transformations between the ECl (Earth Centred Inertial system) 
and the ECEF (Earth Centred Earth Fixed system). The orbital parameters of the satellite con- 
stellation are expressed in the ECEF. To realize the ECEF a network of ground stations receiving 
the GNNS signals has been installed. The GPS uses the World Geodetic System 1984 (WGS-84) 
and Galileo will use the Galileo Terrestrial Reference Frame (GTRF) ( , Altamimi|2009[ l . These 
global reference frames are fixed to the Earth (via ground stations) so their precision and stability 
in time are limited by our knowledge of the Earth dynamics. The main effects are plate tectonic 
motions, tidal effects on the Earth's crust and variations of the Earth rotation rate. In the WGS-84 
the best accuracy achieved is 30 cm ( NIMA|2000| l, with an average stabihty of 4 cm/year ( Al 



|tamimi|2009 ). The use of other space geodetic techniques - VLBI, SLR and DORIS - is necessary 
to achieve a high precision ECEF. The International Terrestrial Reference Frame (ITRF), main- 
tained by the International Earth Rotation and Reference Systems Service, combine efficiently 
these four techniques to reach a stability of about 1 mm/year, inferior by a factor of 10 to science 
requirements in geology ( Altamimi,2009|). 

These considerations led Bartolom CoU ( 2003 | l to propose the project "Systeme de Position- 
nement Relativiste" (SYPOR), i.e. Relativistic Positioning System, an alternative to the scheme 
of usual positioning systems. The idea is to give to the constellation of satellites the possibil- 
ity to constitute by itself a primary and autonomous positioning system, without any a priori 
realization of a terrestrial reference frame. The relativistic positioning system is defined with 
the introduction of null coordinates. They have been reintroduced recently by several articles 
( |Coll and Morales] 199 l[|Rovelli|2002l|Blagojevic et al.|2002[|Lachieze-Rey|2006p . They have 
different names in the literature: "null coordinates", "emission coordinates", "GPS coordinates", 
"GNSS coordinates". In this article we use the first name which is a reference to their geometrical 
properties (see Coll and Pozo (20061 for a detailed article). 

The definition of these coordinates is rather simple, but they are a very powerful tool in 
general relativity. Let us define four particles a = 1,2,3,4 coupled to general relativity. Their 
worldlines Ca are parameterized by their proper time t°. We choose a random origin for t" = 
on each worldline. Let P be an arbitrary event. Then the past null cone of P crosses each of 
the four worldlines in t° (see Figjlji. The quadruplet (t^ ,t^,t^ ,t*)p constitutes the complete 
set of null coordinates of the event P. The protocol to define null coordinates can be seen in a 
different way. The worldline Ca of the particle a defines a one-parameter family of future null 
cones, which can be parameterized by proper time ta (see Fig|2]). The intersection of four future 
null cones t" from four worldlines Ca defines an event with coordinates {t\t^,t^ ,t'^). A user 
receiving these signals knows its position in this particular coordinate system. 

Coll and collaborators ( |Coll et al. ,2006b a ) studied relativistic positioning systems in the case 
of a two-dimensional space-time for geodesic emitters in a Minkowski plane and for static emit- 
ters in the Schwarzschild plane. A relativistic positioning system has been studied in the vicinity 
of the Earth: calculations were performed to first order in a Schwarzschild space-time (Bahder| 
2001 Ruggiero and Tartagha||20(j8 1. |Bini et al.|P008|l described null coordinates in the local 



Fermi frame of a user A "galactic reference system" has been studied, where timing signals 
received by four pulsars were considered as null coordinates ( CoU and Tarantola,2003,,Tartaglia| 




Figure 1 : Defining null coordinates with the null past 
cone of a space-time event: let P be an event in space- 
time. Ca is the worldline of a test particle a parameter- 
ized by its proper time t° ; its origin O is in t" = 0. The 
past null cone of the event P cross Ca at the proper time 
T° . With four different particles with the worldlines Co 
(a = 1, 2, 3, 4), the past null cone of P crosses the four 
worldlines in rj,, rj,, tJ, and t^. Then (r' , r^, t^, t'')p 
are the null coordinates of the event P. 




Figure 2: Defining null coordinates with four one- 
parameter families of null future cones: the worldline 
Ca defines a one-parameter family of null future cone. 
These null cones are hypersurfaces of r° = constant. 
The intersection of four null cones defines the event 
with coordinates (T',r^,r^,r*). The four particles can 
be chosen as four satellites of a constellation of satel- 
lites, broadcasting their proper time. 



et al.||2010[l. The next generation of GNSS could have cr oss-link capabilities (Directorate of 



the Galileo Programme and Navigation Related Activities||2009"| ). Each satellite will broadcast 



proper time to the other satellites in view, as well as their proper time. With this information, one 
could in principle map the space-time geometry in the vicinity of the constellation of satellites 
by solving an inverse problem ( jTarantola et al.|20"09 1 



In this article, we assess the practical feasibility of a relativistic positioning scheme for a 
constellation of satellites in orbit around the Earth. The Schwarzschild geometry is chosen as a 
first approximation of the space-time geometry around the Earth. In section |2] we give solutions 
of geodesic equations for closed time-like orbits - satellites orbits - and for scattering light-like 
orbits - photons orbits. In section |3] we describe two algorithms: (i) how to deduce null coor- 
dinates of an event in space-time knowing its Schwarzschild coordinates and orbital parameters 
of a constellation of four satellites and (ii) how to calculate Schwarzschild coordinates of an 
event knowing four proper times sent by four satellites, and their orbital parameters (reverse al- 
gorithm). We study the speed and accuracy of these two algorithms implemented in Fortran 95. 
Finally, in section]?] the effects of non-gravitational perturbations on the relativistic positioning 
system are assessed. 
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2. Geodesies in the Schwarzschild space-time 



The Schwarzschild space-time is a good approximation of the geometry around the Earth. 
The local inertial coordinate system is tied to the center of the Earth and oriented into 4 mutually 
orthogonal directions f, X, Y, and The Schwarzschild space-time is usually represented by 
the metric in spherical coordinates t, r, 0, and (p, such that X = r sin cos (p, Y - rsinO sin 0, and 
Z - r cos . In these coordinates the metric is 



d/ 



-C(r)dt^ + C-\r)dr^ + r^dQ^ 



where 



C(r) 



2M 
r 



(1) 



(2) 



and M is the Earth mass. We use natural unit^c - G = 1. Geodesies 
are governed by the Hamiltonian: 



H 



-C-\r)p^ + C(r)pj + 



pI + 



(3) 



which admits 8 constants of motion ( Cadez and Gomboc|l996 1: value of Hamiltonian (H) and 
Lagrangian (L), controlUng the relation between time and distance, energy (E - p,), three com- 
ponents of angular momentum (/), longitude of periapsis (w) and time of periapsis passage (tp). 

It is convenient to introduce another local inertial (right-handed) orthonormal tetrad «, and 
62 (c.f Fig.|3]l. « is the constant unit vector pointing in the direction of the angular momentum: 
I = In. The two unit vectors ej and 62 in the orbital plane are such that points in the direction of 
the initial perigee. The components of these vectors with respect to the local Cartesian coordinate 
basis are expressed as: 



ei - (cos CO cos f2 - cos i sin a> sin f2, cos <y sin O + cos t sin lo cos O, sin i sin co) 
€2 - (sin (JL> cos Q. - cos l cos o) sin Q., - sin w sin Q + cos t cos u cos Q, sin l cos oj) 
h - (sin L sin Q, - sin l cos Q., cos l) , 



(4) 



where Q is the longitude of the ascending node and l is the incUnation of the orbit with respect 
to the X - Y plane. 

The only p arameter changing along th e orbit is the true anomaly A. It obeys the differential 



orbit equation ( Cadez and Gomboc 



1996 1: 



du 
dA 



: VA2-m2(1 -u) + B(1-u), 



(5) 



where u - 2M/r, and A - 2ME/1 and B - 2H(2M/l)^ are two constants of motion related to 
orbital energy and orbital angular momentum. After (|5]l is solved for m as a function of A, the 
orbit can be described with 

2M 

r(A) - (ei cos A + €2 sin A). (6) 



'The Earth's local inertial coordinate frame is precessing with respect to the global inertial frame tied to distant stars. 
The precessions are due to different gravitational perturbations of Earth's multipole moments, gravitational perturbations 
of the Moon, Sun and planets. These motions can be well modelled, but are not the matter of this article. 

^To recover usual units one can replace M = Gm^/c^ when measuring distance and M = Gm^/c^ when measuiing 
time, where is the mass of the Earth in the usual units. 

4 




Figure 3: The orbital plane in equatorial coordinates: it unit normal, i inclination, f2 longitude of the ascending node, oj 
longitude of periapsis and A true anomaly. 



The spherical coordinates and if> along the orbit are expressed a^ ( [Cadez and Gomboc 
[T9961 Eqs. 6 - 16): 



tan 



cos 6 - sin L sm(A + w) 
(p - n cos L sm(A + a>) 



2 sin 6* + cos(/l + a>) 
Time and proper time obey the following differential equations: 
df _ 2MA 



dr _ 2MA 
du 



1 



E u2 



u^^JA^-uHl -u) + B(l -u) 



(7) 
(8) 

(9a) 
(9b) 



The differential equations (|5]l and ( |9al ) are formally the same for light-like (where B - 0) 
and time-like orbits. However, solutions depend on the type of orbit, e.g. closed, scattering or 
plunging. In the following subsections we give solutions for closed time-like orbits and scattering 
light-like orbits, which can be used to model satellites and photons trajectory in a GNSS. 

2.1. Time-like geodesies 

GNSS satellites are on closed time-like orbits. After solving (jSj) for such case, we obtain the 
orbit equation 

u{X) = t/2 - (t/2 - I K(m«) + ^ 



2n„ 



(10) 



^Note that sin e = + Vl - cos^ 6». 
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Figure 4: Left: A time-like orbit witli semi-major axis a = 500 Tg and eccentricity e = 0.3. Middle: Time (black) and 
proper time (red) for the same orbit. Right: The difference between time and proper time for the same orbit. All values 
of r, t, and r are in units of M. 



where U\, U2 and U-}, are the roots of the polynomial P(u) - - u^(\ - u) + B{1 - u), and Ua 
and nia are functions of them. K and cn are the complete elliptic integral of the first kind and 
Jacobian elliptic function respectively. The true anomaly A is defined as in the Keplerian case; at 
periapsis it has values Ap = 4n„K{ma)k, where k is an integer. U2 and are related to the radii 
of the periapsis rp = 2M/U2 and apoapsis - 2M/Ui respectively. For circular orbits U2 - t/3, 
and the orbit equation reduces to r = const. - r^- rp. 

Schwarzschild time and proper time alo ng the orbit are obt ained by integrating equations 



- (Nbb. Their expressions can be found in Cadez et al. (2010 1 



A solution of the time-like geodesic equation is illustrated in Fig.|4] Orbital parameters were 
intentionally chosen such that relativistic effects are clearly visible, i.e. periapsis precession and 
the difference between time and proper time. For the values of orbital parameters of the Galileo 
satellites, the difference between coordinate time and proper time goes up to 10 i-is per orbit. 

2.2. Light-like geodesies and ray-tracing 

In the eikonal approximation, electromagnetic signals sent by satellites follow null scattering 
geodesies. After solving (jSjl with B = for such cases, we obtain the orbit equation 

A 



u(A) - U2 - (u2 - uj)cn^ ^K(m) H — nij , (11) 



where the true anomaly A takes values on the interval A e iF(xmax\tn) - K(m), F(^„„„|m) - K(m)). 
Here Xmm = arccos ( Vm2/(m2 - M3)), Xma.^ - arccos Vm2/(m2 - M3)), F is the elliptic integral of 
the first kind, and ui, U2, M3 are the roots of the polynomial P(u) - - u^{\ - u). The constants 
M], U2, M3, m, n depend only on one constant of motion - A - 2ME/1, which is the inverse of the 
impact parameter. 

Coordinate time and proper time as a function of A is obtained by integrating equations (|9ajl 



and (9b 1 (Cadez et al. 2010). In the case of GNSS, the gravitational field is very weak so photon 



orbits are essentially straight lines. Differences between the relativistic and non-relativistic time- 
of-flight are of the order of 1 ns when considering a signal travelling from a satellite to the user 
Calculating null coordinates of an event in space-time requires one to be able to follow light 
rays emitted from the satellites towards that event. Thus, one would like to find the quickest 
way to determine all the constants of motion of a null geodesic that connect a given initial point 
Pi = (f,, r,, 6i, ipi) and a final point !P/ = (f/, rj. Of, iff). ^Cadez and Kostic, ( |2005| ) have proposed a 
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Figure 5: Determining a light-like geodesic between points Pi and ff from known A/1, r,, and rf. The calculated orbit 
is shown in red. This case corresponds to a strong gravitational field in order to see the effect. 




Figure 6: Determining null coordinates from Schwarzschild coordinates. The signal is received at The point of 
emission P, has to be determined by connecting the two points with a light-like geodesic (diagonal line). For clarity, only 
one spatial coordinate is shown. 



ray-tracing method which completely solves this problem in the Schwarzschild space-time. The 
method uses spherical trigonometry and the analytic expression of the orbit equation (111. An 
example of orbit determination from A/l, r,, and r/ is shown in Fig.|5] 



3. Simulation of the constellation of satellites 

3.1. Determining null coordinates 

We assume that the constants of motion of all satellites are known. Their trajectory can 
be calculated and parameterized by their true anomaly X. Let us define a user at the event 
Vo - (k],^o,yo,Zo). He receives signals from four satellites that were sent at !P/ = {ti,Xi,yi,Zi), 
corresponding to /I,- (/ = 1 , 4). Null coordinates of the user at 'Pg are the proper times r,(/l,) of 
the satellites at !P, (see Fig.|6|. We calculate /I, at the emission point !P, using the equation that 
connects To and !P, with a light-like geodesic 

t„-ti(Ad^Tj(ii,(A,)Jo), (12) 
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where ^, - (X,, Yi,Zi) and ^„ - (Xg, Yo,Zo) are respectively the spatial vectors of the satellites 
and the use r. The function Tf calcu lates the time-of-flight of photons between and !P,- as 



described in Cadez and Kostic (2005 1. This equation can be very efficiently solved with known 
algorithms, e.g. Newton's method. Once the value of /i, is determined, it is straightforward to 
calculate the proper time of emission t, for each satellite and thus obtain null coordinates of the 
user at !Po = (ti,T2,T3,T4). 

3.2. Calculating Schwarzschild coordinates 

Here we solve the inverse problem: calculate Schwarzschild coordinates of the event Pg 
from (ti,T2,T3,T4) sent by the four satellites. As we assume that the constants of motion of 
all satellites are known, we can deduce their space-time positions Vi - (f,, jc,, y,, z,) from proper 
times T;. Events !P/ = (?,■, x;, y;, zi) and V,, - (to, Xo, yo, Zo) are connected with light-like geodesies. 
In a flat space-time, they solve the four equations 



to - ti = ^|{.Xi - x„)2 + (j. - yjz + ^j^. _ x^)2_ (13) 

These four equations can be solved for ito,Xo,yo,Zo) by a geometric construction. Let Rj - 
(Xi, Yi,Zj) be the spatial coordinates vectors of the satellites at Vi- The situation is illustrated in 
Fig. [7] where the four green points represent the four satellites at ^, and the red point is the user 

The 4 spheres centred at have radii (?„ - ?,). Thus, the user is at the intersection of the 
four spheres. To find his position, we proceed as follows: suppose that the user's dead reckoning 
coordinates are {t''o\xo\y^o\:^o^). The radii of the four spheres centred at Ri would then be 
(tf^ - ti). In general, these four spheres have no common point. However, the dead reckoning 
position can always be chosen in such a way that any two spheres intersect. Consider the planes 
defined by the circle of intersection of sphere 1 and 2 and sphere 3 and 4. These two planes 
generally intersect along a straight line; call it line 1 . Next consider the intersection of spheres 1 
and 3, and spheres 2 and 4. The corresponding planes intersect along a straight line called line 



{tf\ Xg' , 3'o"\ Zo'O y^*- solution of (|13b, the lines 1 and line 2 generally bypass each other 
We calculate the positions on both lines where they meet at closest distance. The geometrical 
center of the two positions is taken as the better approximation for the spatial position of the 
user R^g \ and the distance between the two points of closest passage (dP) is a measure of how 
close we are to the solution. In the next step, we repeat the procedure with a slightly diff'erent to 
to find the derivative of dP with respect to to. This derivative is then used in Newton's method 
to find a new fj,'' with a smaller dP. With this value a new spatial position P^,^' is calculated 
and the procedure is repeated until dP becomes less than the value prescribed by the accu racy 



2. If ( 13 1 is satisfied, then line 1 and line 2 intersect at the position of the user. However, since 



,(0) ,(0)^ 



13i 



requirement. After n steps (n ~ 5), the procedure converges to {fo"',P["*), the solution of 
Note, however, that a unique solution exists only if the normals to all planes whose intersections 
define lines 1 and 2 do not lie in a plane, i.e. if the four satellites do not lie in one plane. 

However, eq. ( [T3| ) does not take into account the gravitational time delay and the bending 
of light, so our result is not exact. For realistic positions of Galileo satellites the position error 
is a few 10 orbital radii, i.e. a few centimetres. The final correction to the position is made 
by solving the relativistic equations of light propagation from the satellites to the user in the 
following form: 

to-ti^Tf(Ri,Ro), (14) 
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Figure 7: The four green points represent the four satellites at Rj and the red point is the user. 

where Tf(^i,^„) calculates the coordinate time for the light to travel from ^, to Taking 
as an initial approximation, we use a generalization of the classical stellar navigation 
solution to solve ( [T4| . This equation is written in the form: 

= TfiRi, + V„Tf{Ri, ■ + O(d^), (15) 

where Vq is the gradient with respect to the position of the user Rg. We now assume that the 
gravitational field is weak, i.e. \M/R\ <K 1. Then 

Tf{Ri,R„) = \R,-Ri\+0(M) 

and 

\Ro-Ri\ 

= Ui+0(M/R), (16) 

where m, is the unit vector pointing from satellite / to the supposed position of the user. Eq. ([T4]l 
now becomes a set of four linear equations for df and the three components of dR: 

f("> - f,- - TfiitiJ^^^) = -df + Cf!"^ dit. (17) 

The error of the corrected position decreases by 6 to 9 orders of magnitude - bringing the Galileo 
position error to the order of micrometers. If necessary, this error can be decreased even further 
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Figure 8: The orbits of 4 satellites around the Earth. The satellites at initial positions are marked with black dots, and the 
user on Earth's surface with a red dot. The green sphere represents the Earth. The sizes of satellites are not to scale. 
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Table 1: Orbital parameters for 4 satellites: longitude of ascending node (Q.), longitude of perigee (ai), inclination (i), 
major semi-axis (a), eccentricity (e), and time of perigee passage (tp). 



by solving ( 17 1 again with fj,"-^'' 



f!,"' + dt and 



R' 



'(n+l) 



7?!"' + dR 



(18) 



3.3. Accuracy and speed of algorithms 

The above algorithms were tested in a simulation of 4 satellites 5 , (/ = 1 , . . . , 4) in orbit around 
the Earth communicating with a static user on the Earth's surfac^ The only input parameters are 
the orbital parameters of the four satellites (table [TJ and the coordinates of the user. The orbits 
for the four satellites are illustrated in Fig. [8] The Schwarzschild coordinates of the user are 

r„ = 1.595 ■ 10'' M , 0„ = 43.97° , 0„ = 14.5°. 

The simulation runs in the following way. At every time-step 

f<"> = f*"-'^ -(- Af , « = 2,3..., 

where Af = 6 ■ 10'^ M and n counts the time-steps of the simulation: 

1. We calculate null coordinates (Tj"\ t""\ t""\ t^"') of the user from his Schwarzschild coor- 
dinates (f*"\ A:o""'',yo""'\zo""''), as described in Sec. [TTJ (For« = 1 we choose (f^'^Xo"', Jo"', 
{0,ro,e,„4>o)). 



''The problem of the atmospheric perturbations is not in the scope of this paper 
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n 


e/ 






e. 


430 


1.93758E-32 


-6.15473E-27 


-7.34498E-26 


-4.86318E-26 


431 


-2.99206E-32 


-1.78472E-26 


-7.64415E-26 


1.84489E-25 


432 


-1.46740E-31 


-4.82243E-26 


-1.13881E-25 


7.31360E-25 


433 


-1.72335E-32 


-1.09471E-26 


-7.03668E-26 


1.27280E-25 


434 


-1.01827E-32 


-7.659 14E-27 


-5.69550E-26 


9.278 16E-26 



Table 2: Relative errors of the coordinates as defined in l |19) - j20) , for the last five time-steps of the simulation, which 
is equivalent to more than one orbit of the satellites. 



2. From previously calculated null coordinates (Tj"\ t'^\ t^"\ t^"^), we calculate Schwarzschild 
coordinates (tf\ xf\yf\ zf^) for every satellite S , at its new position t*"' by numerically 
solving the equation T(Af ^) - rf' for Af \ 

3. From these Schwarzschild coordinates (tf \x^"\y^"\z^"^), we calculate Schwarzschild co- 
ordinates (/g\x^o\y^"\zo'^) of the user as described in Sec. 13. 21 

The numerical code of this simulation is written in Fortran 90 and is publicly available on 
the website atlas . estec . esa. int/ariadnet This code can be easily generalized to include 
a moving user, more satellites or communications between all the satellites. 

Accuracy. The accuracy of these algorithms is tested by comparing the initial Schwarzschild 
coordinates of the user with his Schwarzschild coordinates calculated at each time step. As the 
user is static, his coordinates (xo,yo,Zo) should be constant during the simulation. 
In table|2]we show the relative errors of all the coordinates, defined as 

(n) - 4"' 

(«) ^ O ^jg^ 

and 

Jl) An) „(1) In) ,(1) An) 



w) _ An) _ An) _ ~" /oan 

- (1) , Cy - (1) , - (1) ■ K^^) 

Xq yo Zo 

The relative error of the coordinate f is ~ 10"^^, and relative errors of the coordinates x, y, and z 
are few orders of magnitude larger, i.e. ~ 10"^^ - 10"^^. 

The algorithm for determining Schwarzschild coordinates from null coordinates (Sec. 13.21 1 
works only if the four satellites are not in the same plane. This may happen during a simulation in 
which case the resulting position has a very large error or is left undetermined. As an example, we 
show in table [3]relative errors where the positions cannot be determined for a given configuration 
of the satellites, which is reflected in a jump in the relative errors from ~ m-^^ to ~ 10"- for 
the coordinate f, and from ~ 10"^^ to ~ 10^^ for the coordinate z- In the case of GNSS satellites 
such a situation should never occur, since there are always more than four satellites visible and 
it should always be possible to choose four that do not lie in the same plane. 

Speed of calculations. The simulation was tested on a PC with an Intel Core2 Quad CPU Proces- 
sor Q6600 at 2.4 GHz and 4GB RAM. The OS was Linux x86_64 with kernel 2.6.28- 18-generic, 
with the Intel Fortran Compiler 10.1. The maximum number of time-steps is n„,ax - 434, equiv- 
alent to more than one orbit of the satellites. The calculation for steps 1 to 3 described above 
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10 


-5.15060E-02 


-7.85972E+02 


-7.40892E+01 


-2.31309E+03 


11 


-1.51029E-02 


-2.5671 lE+02 


-2.44827E+01 


-7.5624 lE-i-02 


12 


-9.12949E-03 


-1.71063E+02 


-1.66154E+01 


-5.04159E-I-02 



repeated 434 times (i.e. for 434 time-steps) takes 67.6 seconds to execute (0.1558 seconds/time- 
step). The time of calculation for steps 2 and 3 (repeated 434 times) is 26.6 seconds (0.0682 
seconds/time-step). These two steps are actually the ones that the users will perform to deter- 
mine their positions. 

The simulation was also tested on a laptop with an Intel Core2 Duo CPU Processor P8600 
at 2.4 GHz and 4GB RAM. The OS is the same as before, but with the Intel Fortran Compiler 
11.1. In this case, some processor-specific compiler optimizations are enabled. The maximum 
number of time-steps is the same as before: n^ax = 434. The time of calculation for steps 1-3 de- 
scribed above repeated 434 times (i.e. for 434 time-steps) is 61.9 seconds (0.1426 seconds/time- 
step). The time of calculation for steps 2 and 3 (repeated 434 times) is 26.5 seconds (0.0588 
seconds/time-step). 

In both cases, a 25 - 32-digit accuracy is achieved when determining the space-time position 
of the user. 

4. Effects of non-gravitational perturbations 

The important non-gravitational perturbations are those governed by stochastic noise, i.e. by 
phenomena that cannot be predicted in advance. Clock noise, solar radiation pressure, solar wind 
pressure and collisions with interplanetary dust particles are sources of such noise. The nature of 
these forces and also the nature of the uncertainty in knowing their magnitude and direction are 
quite difl'erent. 

Clock errors. Clock errors affect the position determination by giving false information on the 
time of flight from the satellite to the user. Assuming that a typical clock error can be described 
as a flicker noise with a standard deviation error of 1 ns/day, we can assign a displacement error 
of dt =0.3 m per day to stochastic clock perturbations. 

Solar radiation. The solar radiation pressure produces a force 



where f is the radius vector from the Sun to the satellite, ^ is the cross section of the satellite 
from the impinging direction, a is the effective albedo of the sateUite (generally depending on 
the angle of incidence of the solar light on the sateUite) and /; is a tensor attached to the sateUite, 

measuring the effective momentum of reflected Ught. Assuming that a, rj and J?l are constants. 



Table 3: Relative errors for a planar configuration of satellites. 




(21) 
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Figure 9: As the Earth (blue sphere) moves about the Sun (yellow sphere), the solar radiation pressure makes the satel- 
lites' orbital angular momentum precess so that the tip of the orbital angular momentum follows the red ellipse; the gray 
circle shows where the tip of angular momentum would be without precession. 



the effect can be calculated very precisely since the solar luminosity is a well known constant. In 
this case the main effect of solar radiation pressure is a precession of the satellite orbital angular 
momentum with an angular velocity 

Lojrs m 

where n is defined in section|2] m is the mass of the satellite, Qyeai is the angular velocity of the 
solar radiation force with respect to the orbital angular velocity of a satellite, i.e. fiyear - f^^' 
a>s = 27T Ps (where Ps ~ \ day) is the orbital angular velocity of the satellite and is the distance 
of the satellite from the center of the Earth. This effect, schematically shown in Fig. |9] makes 
the orbital angular momentum oscillate with the orbital period of the Earth around the Sun. The 
displacement of the satellite's position due to solar radiation pressure as a function of time can 
be written as 

'Irp = y-y-P^s sin (fiy^rf) cos(w,f + 6). (23) 

The amplitude of this term can be considered as a measure of the solar radiation strength. If 
we assume the satellite to have a mass of 600 kg, a cross sectional area of 0.5 m^ and take 
- 1300 W/m^ for the solar constant, we obtain drp ~ 0.17 m. 

A realistic calculation of the effect of the solar radiation pressure is somewhat more compli- 
cated, since a, rj and ^ depend on the orientation of the satellite with respect to the direction to 

the Sun. It is probably safe to assume that this dependence can be known to better than 1 % or 
maybe even 0. 1 % in which case the unpredictable stochastic part of the effect of solar radiation 
pressure would be always below millimeters. 

Solar wind. The perturbing force due to solar wind pressure is similar to the radiation pressure 
force since it originates from the same direction. However the pressure term ^l^^'must be 

replaced by 4^v, where M and v are respectively the solar wind mass loss rate and velocity. 
Their typical values found in literature are: M ~ 4 x lO'** g/year and v ~ 500 km.s"'. Using 
these numbers we deduce that solar wind pressure is about 200 times weaker than solar radiation 
pressure. 
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Because of this weakness it may not be important at this point to be able to determine the 
uncertainty in the solar wind albedo and the equivalent of 77. It is certainly not as important as 
finding the truly stochastic part of solar wind, which is constituted by the coronal mass ejections 
hitting the satellite. Coronal mass ejections are violent outbursts from the solar surface, occ urring 



every few days and carrying away masses 10 to 10 g at speeds up to 5000 km.s (Landi 



et al.pOlQ) Zhang et al. 2010| l. We assume that the coronal mass flow is spread evenly in a 



spherical shell moving away from the Surj^ When particles from this shell hit the satellite, it 
gains momentum 6p. The displacement after one day is 

5p 

d,^e = — ■ 1 day (24) 
m 

Using the largest numbers quoted, we obtain an estimate dcme ~ 0.0003 m for the purely stochas- 
tic part of solar wind displacement a day after the coronal mass ejection. 

Interplanetary dust. A relatively small non-gravitational perturbation, but completely stochastic 
in nature is produced by collisions of satellites with interplanetary dust. Our estimates of these 



perturbations are based on data from Nesvorny et al. ( 2006 ), who studied the dust bands Karin 



and Veritas. The authors claim that these dust bands contribute from 30 to 50% of all inter- 
planetary dust in the near Solar system. They contribute 15000 to 20000 tons per year in dust 
accretion rate to the whole Earth. We assume that in the vicinity of the Earth the dust accretion 
is approximately isotropic and that it is proportional to the accreting area {AnR^ - AJ[ - 2m^). 
Taking 36000 tons per year as the total mass accretion rate on the Earth, we thus estimate the 
dust accretion rate per satellite to 1.4 x lO"'* g/year. 

We come to the following main conclusions: solar system dust presents a very mild drag 
resistance to satellite motion around the Earth. The typical orbital decay time scale is of the 
order of ~ 3 ■ lO^years. The stochastic component is contributed mostly by collisions with 100- 
200 pm dust particles which move with respect to the Earth (and with respect to satellites) with 
an average velocity ~ 17 km.s ' and occur with a probability of less than one per year. The 
velocity change due to such a collision is Sv^ust ~ 5 - 20 ■ 10"** m.s"'. The displacement of a 
satellite due to such a collision observed after one day would amount to t/dujr - 5 - 20 mm. Note 
that such collisions could be detected and their impact measured if the satellites would monitor 
their mutual positions, since the probability for more than one collision to occur during the same 
day is very small. 

Comparison. It is difficult to compare the strength of these so different non-gravitational pertur- 
bations without more extensive data on their character, changing strength and without carefully 
considering statistical properties of these noise sources. Such a study would go beyond the scope 
of this article and may not even be very useful before it could be verified experimentally. To 
make a very rough comparison of importance of each of the listed perturbations, we assign to 
each perturbation a "force" and a displacement/day. As the worst case we simply use the full 
perturbing force to calculate the displacement of a satellite after one day. This is a gross overes- 
timate in the case of solar radiation pressure or solar wind pressure. We attempt to refine such 
noise upper limit by estimating only the stochastic part of 1-day position error, i.e. the error 
remaining after the known part of the perturbing force has been taken into account with the best 



This is certainly not a realistic assumption with respect to each single event, but may not be so bad on long term 
average. 
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Table 4: Non-gravitational noise sources 



perturbation 


magnitude 


displacement/day 


stochastic component 


clock noise 


Ins/day 




0.3m 


radiation pressure 


4.3xlO-^Nm-2 


13m 


<0.001m? 


solar wind pressure 


2.2xlO-'*Nm-2 


0.007m 


~ 0.0003m? 


dust collisions 


~ lOOyum part 




~ 0.00005m 



available data. Dust collision differs from all other perturbations considering their very low prob- 
ability rate. Thus their displacement per day is just a very long term average. Our estimates are 
given in Table[4] The reader should be aware of the large uncertainty involved in these estimates. 

Even if the above estimates are very rough it is clear that timing noise is by far the most 
important noise source that must be controlled on a longer time scale. However, it is worth 
pointing out that the non-gravitational perturbations, apart from the slow drag mentioned above, 
do not systematically change the orbital momentum of the satellites (after all the precessions due 
to gravitational perturbations have been properly taken into account). It implies that the area 
swept by the radius vector from the center of the Earth to each satellite can be a clock with a 
very good long term phase stability, limited in principle only by the uncertainty in the swept 
area. If each satellite was capable of detecting the timing signal of all the other visible satellites, 
it would be possible to use the timing capability of the constellation of satellites to determine its 
own space-time position. 

5. Conclusion 

We studied the use of null coordinates. They are user independent and realized by radio wave 
communications from a constellation of four satellites to a GNSS user. We introduced a local 
Schwarzschild coordinate system with three spatial orthonormal directions and one time-oriented 
vector, with coordinates X, Y, Z, t. Before taking into account the gravitational perturbations due 
to the Moon, Sun, planets, obliquity and rotation of the Earth, etc. . . , this local coordinate system 
can be considered as a global inertial system, i.e. oriented in fixed directions with respect to dis- 
tant stars. Since the main effect of gravitational perturbations is to make the local Schwarzschild 
frame precess about different axes determined by the orbits of the perturbing bodies, and these 
precessions are very slow, the local Schwarzschild coordinate frame can to a high degree of ac- 
curacy be considered as inertial. Due to this fact it is possible to decouple the problem of the 
motion in the local Schwarzschild frame from the problem of connecting the local Schwarzschild 
frame to the global inertial frame. This second problem is well understood in the framework of 
classical non-re lativistic gravitational perturbation theory, and will be the subject of a following 
article. 

In this article we concentrated on finding algorithms to describe the dynamics of bodies in 
the local Schwarzschild frame in full general relativity, in defining null coordinates in space-time 
that are tied to the constellation of satellites and in reading these coordinates in order to deter- 
mine the Schwarzschild coordinates of an event in space-time. We obtained the analytic solutions 
for light-like (scattering) and time-like (closed) geodesies. We defined two algorithms: (i) an al- 
gorithm that calculates null coordinates (t\,T2, T3, T4) corresponding to the local Schwarzschild 
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coordinates Xo,yo,Zo) of a user, and (ii) the "reverse" algorithm that calculates space-time co- 
ordinates (to, Xg,yo, Zo) of a user from its null-coordinates (tj , T2, r3, T4). We assume that orbital 
parameters of satellites are known. 

The two algorithms have been written in Fortran 90 and checked against each other for con- 
sistency. All codes have been optimized for weak gravitational fields, and tested on a PC with 
a quad-core CPU and a laptop with a dual-core CPU. It takes ~ 60 - 70 ms to determine all 
four coordinates with 25 - 32-digit accuracy. This makes our codes feasible and could be used 
in modem positioning devices. The effects of non-gravitational perturbations have been studied. 
We focused mainly on the stochastic part of these perturbations. Clock noise has been identified 
as the most important contributor affecting the accuracy of position determination. By magnitude 
the second perturbing force is solar radiation pressure. We note, however, that this force is very 
stable and could, at least in principle, be modeled to a very high precision. Therefore, its stochas- 
tic effect is not expected to contribute significantly more to position noise than the remaining two 
truly stochastic perturbations: solar wind pressure and collisions with interplanetary dust. 

We believe to have shown that the use of a fully relativistic code in GNSS systems offers 
an interesting alternative to the use of post-Newtonian approximations, and presents no technical 
obstacle. The advantages of using a fully relativistic code over the use of classical post-newtonian 
codes go well beyond aesthetics. The stochastic part of the non-gravitational perturbation forces 
has no preferred direction; therefore they are not expected to change average values of orbital 
angular momenta of satellites. It implies that area velocities of satellites are good constants of 
motion perturbed only by deterministic gravitational perturbations. Therefore we suggest that the 
area covered by the radius vector from the center of the Earth to the satellites be used as a stable 
measure of proper time, in particular to correct the long term clock drifts. The constellation 
would have a greater autonomy in defining its own local Schwarzschild frame in which the 
predictability of dynamics can be taken advantage of 

This concept of autonomy offers an increased accuracy and long term stability to the global 
positioning system. It promises us a tool for a systematic study of the space-time metric around 
the Earth and a deeper insight into the Earth internal dynamics. In order to realize such an au- 
tonomous constellation, it must be self-consistent. By this we mean that if satellites determine 
their orbits by using positioning data from other satellites in the constellation, they should obtain 
the same constants of motion that were used in the definition of the global positioning system. 
Such self-consistency at the maximum level of precision is not obvious, but can be reached if 
each satellite is also a receiver of timing signals of all other satellites. We believe we can de- 
sign a mathematical procedure to drive any initial GNSS solution toward the best self consistent 
solution. 
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